
# library(pacman)
# p_load(dplyr, ggplot2, countrycode, igraph, Matrix, tidygraph, ggraph, DBI, RSQLite) 
library(dplyr)
library(ggplot2)
library(countrycode)
library(igraph)
library(Matrix)
library(tidygraph)
library(ggraph)
library(DBI)
library(RSQLite)

con <- dbConnect(RSQLite::SQLite(), "iso_standards.sqlite") # Attaching to database

participation <- dbReadTable(con, "participants") # Extracting TC-membership dataset

dbDisconnect(con) # Disconnecting from database

get_network <- function(year){ # Function to create a network for each year of observation
  
  participation <- participation %>%
    filter(year == !!year) # Select year from participation dataset
  
  edgelist <- participation %>% 
    filter(membership %in% c("P-member", "Secretariat")) %>% # Select only countries that were P-member or in Secretariat
    select(country, title) %>% 
    na.omit()
  
  graph <- graph_from_data_frame(edgelist, directed = FALSE, vertices = NULL) # Make a graph from the dataframe 
  
  # Making the graph bimodal (i.e. from one country to another, not from one country to a TC): Countries become nodes, TCs become edges
  types  <- bipartite.mapping(graph)$type # Find bipartite node types
  matrix <- as_incidence_matrix(graph, types = types) # Make matrix out of dataframe
  mode(matrix) <- "numeric" # Ensure it's numeric
  diag(matrix) <- 0 # 0 at the diagonal
  matrix <- as(matrix, "sparseMatrix") # Sparse to make it easier to compute with
  
  matrix  <- tcrossprod(matrix) # Transpose matrix
  
  g <- graph_from_adjacency_matrix(matrix, mode = "undirected", weighted = TRUE, diag = FALSE) # Make from adjacency matrix to edgelist
  
  degree <- degree(g) # Get nodes
  nr <- attributes(degree) %>% as_tibble() # Add variable attributes to the nodes (country characteristics, e.g. which continent they belong to)
  degree <- degree %>% as_tibble() %>% # Make into a dataframe and calculate degree and eigenvector
    mutate(member = nr$names) %>%
    mutate(degree_total = sum(value, na.rm = TRUE)) %>%
    mutate(degree_percent = value/degree_total*100) %>%
    mutate(year = as.numeric(year))
  
  eigen_centrality1<- eigen_centrality(g)
  eigen_centrality <- eigen_centrality1[[1]]
  country <- attributes(eigen_centrality) %>% as_tibble()
  
  eigen_centrality <- eigen_centrality %>% 
    as_tibble() %>%
    mutate(member = country$names) %>%
    mutate(eigen_centrality_total = sum(value, na.rm = TRUE)) %>%
    mutate(eigen_centrality_percent = round(value/eigen_centrality_total*100, 3),
           value = round(value, 3)) %>%
    mutate(year = as.numeric(year))
  
  return(list(g, degree, eigen_centrality)) # Return the graph, the dataframe with degrees and the dataframe with eigenvector
  
}

gdata <- get_network(2022) # Make the network for 2022

continents <- participation %>% # Make dataframe with each continent that countries belong to 
  mutate(Continent = countrycode(sourcevar = country,
                                 origin = "country.name",
                                 destination = "continent")) %>%
  mutate(Continent = ifelse(country == "Serbia and Montenegro", "Europe", Continent)) %>% # Had to be added manually
  select(country, Continent) %>% 
  unique()

gdata[[1]] %>%  # The first element in the function spits out the graph
  as_tbl_graph() %>% # Turn graph into tidy format and plot
  activate(nodes) %>% # Working with the nodes
  left_join(continents, by = c("name" = "country")) %>% # Add continents to data
  mutate(name = ifelse(name == "United Kingdom", "UK",
                       ifelse(name == "United States", "USA",
                              ifelse(name == "Russian Federation", "Russia", name)))) %>% # Shorter names to better fit into plot
  mutate(degree = gdata[[2]]$value) %>% # Add degree to show how many committees each country participates in
  mutate(`Eigenvector centrality` = gdata[[3]]$value) %>% # Add eigenvector to show how influential a country is in the network
  ggraph(layout = "stress") + # Type of graph layout (looks nicest)
  geom_edge_link(alpha = 0.01) +
  ## REPLACE WITH LINE 79-83 IF RUNNING ON A STRONG COMPUTER TO REPLICATE EXACTLY
  # geom_edge_link(
  #   alpha = 0.1,
  #   aes(edge_width = weight,
  #      color = ifelse(weight >= 80, "high", "low")), # Add color code to the edges to show how central nodes have more connections
  #  show.legend = FALSE) +
  scale_edge_color_manual(values = c("high" = "darkblue", "low" = "grey66")) + # Add color to the color code
  geom_node_point(aes(shape = Continent, 
                      size = `Eigenvector centrality`, # The size of the node becomes the eigenvector. 
                      color = Continent),
                  alpha = 0.8) + 
  # Some specifics to make the plot fit nicely on the page (reduced sizes for the repex, because what looks nice in R doesn't always look nice on the page, and vice versa)
  scale_shape_manual(values = c(0, 1, 2, 5, 6)) + 
  scale_color_manual(values = c("#0F0F0F", "#0F0F0F", "#0F0F0F", "#0F0F0F", "#0F0F0F")) +
  geom_node_text(aes(filter = degree >= 120, label = name), size = 4) +
  guides(color = guide_legend(override.aes = list(size = 10))) +
  guides(size = "none") +
  scale_edge_width_continuous(range = c(0.01, 1)) +
  scale_size_continuous(range = c(4, 8)) +
  theme_graph() +
  theme(legend.position = "bottom",
        text = element_text(size = 20),
        legend.key.size = unit(1, 'cm'))
